from scipy import linalg
import numpy as np
'''
# det函数计算行列式
A = np.array([[1,2,-1],[3,4,7],[9,10,-2]])
x= linalg.det(A)
print(x)

# 计算线性方程组的特征值和特征向量
A = np.array([[1,2,-1],[5,3,4],[9,7,11]])
l,v= linalg.eig(A)
print ('特征值是:',l)
print ('特征向量是:',v)

A = np.array([[1,2,-1],[5,3,4],[9,7,11]])
i= linalg.inv(A)
print ('逆矩阵是:',i)
ii= linalg.inv(i)
print ('逆矩阵是:',ii)

# 矩阵分解
A= np.array([[1,2,-1],[5,3,4],[9,7,11]])
result= linalg.lu(A)
print('矩阵A分解的三个矩阵为:\n')
for i in result:
    print(i)
    print("\n")
r1=np.matmul(result[0],result[1])
rA=np.matmul(r1,result[2])
print('result[0]*result[1]*result[2]=\n')
print(rA)

A = np.array([[1,2,-1],[5,3,4],[9,7,11]])
result= linalg.qr(A)
print('矩阵A分解的两个矩阵为:\n')
for i in result:
    print(i)
    print('\n')
rA=np.matmul(result[0],result[1])
print('result[0]*result[1]=\n')
print(rA)
print('矩阵q的行列式等于:',linalg.det(result[0]))

A = np.array([[1,2,-1],[5,3,4],[9,7,11]])
b=np.array([-1,2,4])
x= linalg.solve(A,b)
print('原方程组的解为x:\t')
print(x)
print('验算 A*x =\t')
print(np.matmul(A,x))

# interpolate 数据插值
from scipy import interpolate
import pylab as pl
x=np.linspace(0,10,11) 
# 创建一个等间隔数值序列
y=np.sin(x)
xnew=np.linspace(0,10,101)
pl.plot(x,y,'ro') # red circle
list1=['linear','nearest','cubic']
# 线性插值，最近邻插值，三次样条插值
for kind in list1:
    print(kind)
    #创建一个一维插值函数
    f=interpolate.interp1d(x,y,kind)
    ynew=f(xnew)
    pl.plot(xnew,ynew,label=kind)
pl.legend(loc='lower right') # 图例
pl.show()

# UnivariateSpline 一维样条插值
from scipy import interpolate
import matplotlib.pyplot as plt
x1=np.linspace(0,10,20)
y1=np.sin(x1)
sx1=np.linspace(0,12,100) # x范围超过[0,10],外插值
func1=interpolate.UnivariateSpline(x1,y1,s=0) # 强制通过所有点
sy1=func1(sx1)
plt.plot(x1,y1,'o')
plt.plot(sx1,sy1)
plt.show()

from scipy import interpolate
import matplotlib.pyplot as plt
x1=np.linspace(0,10,20)
y1=np.sin(x1)
sx1=np.linspace(0,12,100) # x范围超过[0,10],外插值
func1=interpolate.UnivariateSpline(x1,y1,s=0.5) # 不强制通过所有点
sy1=func1(sx1)
plt.plot(x1,y1,'o')
plt.plot(sx1,sy1)
plt.show()

# RectBivariateSpline 二维样条插值
from scipy import interpolate
import matplotlib.pyplot as plt
x=np.arange(-5.01,5.01,0.25)
y=np.arange(-5.01,5.01,0.25)
xx,yy=np.meshgrid(x,y)
z=np.sin(xx**2+yy**2)
f = interpolate.RectBivariateSpline(x, y, z)
xnew = np.arange(-5.01,5.01,1e-2)
ynew = np.arange(-5.01,5.01,1e-2)
znew = f(xnew.ravel(), ynew.ravel())
#插值函数的自变量是网格矩阵，返回的函数f自变量要求是一维数组。
plt.plot(x, z[0, :],'ro-',xnew, znew[0, :],'b-')
plt.show()

# 三维样条插值
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

x= np.arange(-1.5,1.5,0.3)
y=np.arange(-1.5,1.5,0.3)
xx,yy=np.meshgrid(x,y)
z=np.sin(xx**2+yy**2)
fig=plt.figure(figsize=(18,9))
ax1=plt.subplot(1,2,1,projection='3d')
surf1=ax1.plot_surface(xx,yy,z)  
# 绘制原始数据构成的第一个子图
f= interpolate.RectBivariateSpline(x,y,z)
# 插值函数的自变量是网格矩阵，返回的函数f自变量要求是一维数组。
# 采取步长更小的网格
xnew =np.arange(-1.5,1.5,1e-2)
ynew =np.arange(-1.5,1.5,1e-2)
znew =f(xnew, ynew)
xx_new,yy_new=np.meshgrid(xnew,ynew)
ax2=plt.subplot(1,2,2,projection='3d');
surf2=ax2.plot_surface(xx_new,yy_new,znew)
# 绘制原始数据构成的第二个子图
plt.show()

# 多项式拉格朗日插值
from scipy.interpolate import lagrange
import numpy as np
x= np.array([0, 1, 2])
y= x**3
poly = lagrange(x, y)
print(type(poly))
print(poly)

# 最小二乘的多项式拟合
import numpy as np
import matplotlib.pyplot as plt
xx=np.random.rand(50)*4*np.pi-2*np.pi
f=lambda x:np.sin(x)+0.5*x
yy=f(xx)
print(xx[:10].round(2))
print(yy[:10],round(2))
My_fit=np.polyfit(xx,yy,5)
ry=np.polyval(My_fit,xx)
plt.plot(xx,yy,'b^',label='f(x)')
plt.plot(xx,ry,'r.',label='polynomial fitting')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()
'''
#非线性最小二乘法拟合
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
x=np.array(range(10))
y =np.array([0,1,2,3,4,5,4,3,2,1])
def gaussian(x,*param):
    return param[0]*np.exp(-np.power(x-param[2], 2.)/(2 *np.power(param[4], 2.)))+param[1]*np.exp(-np.power(x-param[3], 2.)/(2 *np.power(param[5], 2.)))
#定义模型及参数
popt,pcov= curve_fit(gaussian,x,y,p0=[3,4,3,6,1,1])
print('模型的系数为:',popt)
print('系数的协方差矩阵为:',pcov)

plt.plot(x,y,'b+:',label='data')
plt.plot(x,gaussian(x,*popt),'ro:',label='fit')
plt.legend()
plt.show()